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We investigate the dynamical and steady-state spin response of the nonequilibrium Anderson 
model to magnetic fields, bias voltage, and temperature using a numerically exact method combin- 
ing a bold-line quantum Monte Carlo technique with the memory function formalism. We obtain 
converged results in a range of previously inaccessible regimes, in particular the crossover to the 
Kondo domain. We provide detailed predictions for novel nonequilibrium phenomena, including 
non-monotonic temperature dependence of observables at high bias voltage and oscillatory quench 
dynamics at high magnetic fields. 



Strongly correlated open quantum systems appear in 
a wide variety of physical situations, including quantum 
dots in semiconductor heterostructures P] [2], molecu- 
lar electronics [H H] and the dynamics of cold atoms 
[5]. These systems consist of a finite, interacting region 
coupled to a continuous set of non-interacting "bath" or 
"lead" states which may be maintained at differing ther- 
modynamic states in such a way that no equilibrium 
state can exist. It is natural to describe open systems 
in terms of quantum impurity models, which have been 
used in the description of magnetic impurities in metals 
[5] , the adsorption of atoms on a surface [7J and as aux- 
iliary problems in the dynamical mean field approxima- 
tion to extended lattice systems |H]. More recently, they 
have also been of interest in the nonequilibrium context 
of mesoscopic transport |5J [TD] and nano-systems coupled 
to broad leads [2J. 

While attempts are being made to connect nonequilib- 
rium physics to equilibrium concepts , the nonequilib- 
rium steady state properties of correlated quantum sys- 
tems continue to present a formidable challenge to our 
theoretical understanding. The main difficulty is that 
a rigorous evaluation of the long-time and steady state 
properties requires an accurate time propagation, start- 
ing from some known initial state and reaching all the 
way to the steady state. When this relaxation occurs 
quickly, a range of powerful semi-analytical |12H14| and 
numerical methods |15ff2"2] are applicable. However, dy- 
namics in strongly correlated systems may exhibit a sep- 
aration of timescales — for example, the spin-relaxation 
dynamics in the Kondo regime of a quantum dot are or- 
ders of magnitude slower than those of the corresponding 
charge relaxation. Existing theoretical methods are un- 
able to resolve these timescales reliably. 

In this Letter we show that a combination of bold- 
line diagrammatic Monte Carlo methods [211 123"] and the 
memory-function approach [22; enables us to significantly 
extend the time regime accessible and can, in some cases, 
access steady state information within the Kondo regime. 



The method is numerically exact and provides unbiased 
error estimates. While the calculations presented here for 
the single impurity Anderson model, a minimal model for 
strong interactions in the presence of baths, the method- 
ology is applicable to any quantum impurity model. 

The Anderson impurity model is defined by the Hamil- 
tonian 

H = H S + H B + V, (1) 

where Hg describes the interacting system (or dot) part, 
Hb the non-interacting bath (or leads) part, and V the 
system-bath coupling part: 

Hs=J2 £ A + ^j^H' ( 2 ) 

i=t4. 

Hb = ^ s lk a\ k a ik , (3) 

V= Ukdi(i\ k +t* k a lk d\. (4) 

Here f an d J. represent electronic spin, the di and d\ 
are fermionic system operators for dot states with energy 
£j, ciik and a\ k are fermionic lead operators with energy 
Eik and the U k are coupling constants, k is an index 
iterating over the lead states. Both the eh~ and the ti k are 
fully defined by the system-lead coupling density Y (e) = 

Refs. [2"4Tl2l)] have shown that the reduced density ma- 
trix a (t) = Ttb {p {t)} (p(t) being the full density ma- 
trix and Tr# {...} denoting a trace over all bath degrees 
of freedom) of any system of the form of Eq. ([I]) exactly 
obeys the Nakajima-Zwanzig-Mori equation 

ih^^-=£ Hs <T{t)+#(t)-l [ dr«(r)<r(*-r). (5) 
at h J 

Here, the Liouvillian superoperator Ch s A = [Hs,A] de- 
notes a commutation with the system Hamiltonian Hg. 
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with the same notation defining Ly and Cr] "& (t) is an 
initial correlation term which vanishes for factorized ini- 
tial conditions p(0) = pr <£> f (0); and p_e is the initial 
bath density matrix, k is known as the memory kernel 
and may be obtained by solving [37] 

k (t) = ih® (t) -$(*)£<? + - / dr$ (i - r) « (r) ,(6) 



where the superoperator $ (t) = Tr^ |/Ve R £Ht / os| 
must in general be obtained from a many body com- 
putation whose expense rapidly increases as t increases. 
Evaluation of $ (t) for t up to a cutoff time i c allows an 
exact evaluation of k (t < t c ). Setting k (t > t c ) = de- 
fines the cutoff approximation, whose convergence may 
be monitored from the dependence of results on t c as t c 
is increased. In the case of the Anderson impurity model, 
Ref. [22] has shown that if one is only interested in eval- 
uating the diagonal elements of the density matrix, all 
the supermatrix elements $>ij. qq < = <&\q) (q'\ of 

<I> having i ^ j or q ^ q' can be set to zero, with the 
remaining elements determined by: 

®ii m = S t0 (tp$ + ^ + g a ^(2) _ ^ 

+ 5 i2 + + S« " <P%) ,(7) 

(t) = 2z3^Tr B {p B 4 m) (t) \q)} , (8) 



(i) 



l(2) 



where 

t lk d t d\d ia \ k and A<: 4) 

The evaluation of the ip q q (t) has previously been per- 
formed with real time path integral Monte Carlo (RT- 
PIMC) methods [15] [22j (28j , revealing that, in the pres- 
ence of strong electronic correlations, the memory ker- 
nel may develop long tails. Near the Kondo regime this 
effect becomes particularly pronounced, making it im- 
possible to converge the cutoff approximation and high- 
lighting the need for methods able to obtain k for longer 
times. Here we show that the problem can to a large 
extent be solved by using the bold expansion for impu- 
rity models [53], a technique related to bold-line meth- 
ods for lattice systems |29H31| . The bold expansion is 
based on a stochastic Monte Carlo sampling of diagram- 
matic corrections to the propagators obtained from an 
infinite partial summation, rather than a sampling of all 
diagrams. The resulting procedure converges at lower 
expansion order and greatly reduces the severity of the 
dynamic sign problem, in practice more than doubling 
the accessible time scales. Even with bold methods, a di- 
rect description of the slow spin dynamics remains out of 
reach — however, the bold method does allow converged 
access to such dynamics within the memory formalism. 

In the noncquilibrium case, diagrams must be evalu- 
ated on the Keldysh contour. This had previously been 
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Figure 1. Self energies and vertex equations used within the 
OCA based bold expansion. Solid lines represent bare prop- 
agators, bold lines are dressed propagators, wavy lines are 
hybridization interactions and shaded regions are vertex func- 
tions. The vertices are defined on the unfolded Keldysh con- 
tour, such that the final time on the contour is marked by the 
central "X" and both edges of the contour stand for the initial 
time. The NCA is obtained by taking only the first diagram 
in the self energy and the first two diagrams in the vertex. 



done with "BoldNCA", using the non-crossing approxi- 
mation (NCA) as the underlying partial summation, but 
in the course of this work it has become necessary to em- 
ploy a "BoldOCA" built on the more precise one-crossing 
approximation (OCA) [231 132| . Fig. [TJ illustrates this in 
diagrammatic terms: the bold- line propagators and ver- 
tex functions (which allow for the summation over hy- 
bridization lines connecting pairs of times on the two dif- 
ferent halves of the Keldysh contour) are obtained from 
the solution of the NCA or OCA equations, and used in 
an expansion which samples diagrams of all crossing or- 
ders. Unbiased error estimates are obtained by jackknife 
analysis on data from multiple, uncorrelated Monte Carlo 
runs (typically 5-10). 

We assume a wide, flat band Ti (E) = Tf (E)+Tf (E) 
with (E) — t r= r=n — 7t ; here e c and 

V I (l_|_ e ^(E-e ) Vl_|_ e -«/(E+s c ) y c 

v are the band cutoff energy and its inverse cutoff width, 
and L and R are respectively left and right lead in- 
dices. We restrict our calculations to the symmetric 
Anderson impurity model in a magnetic field h, setting 
Si = — % ± \ (the formalism is more general and does 
not rely on this symmetry). We choose T as our energy 
unit, and throughout the rest of this paper set U = 5T, 
e c = 10r and Tv = 10. The initial conditions are deter- 
mined by assuming an initially decoupled system, having 
left and right leads thermally equilibrated at a tempera- 
ture (3 and chemical potentials p.^ = \ and (1r = 
respectively. This defines the lesser and greater hy- 



\ R ) M = -ifh(R) M T L[R) (w) 



bridization functions A] 

and A> (u>) = i (l — Jl(r) (w)) ^l(r) ( w )j which depend 
on the temperature and chemical potentials through the 

r 2 v- At 

these parameters, the Kondo temperature is given by 

£~ 3.4 |33]. 

In equilibrium, the magnetization can be evaluated ex- 



Fermi occupation function /l(r) (w) 



r/3 
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actly on the Matsubara axis [231 l3~il [55] . As this mag- 
netization corresponds to the steady state magnetization 
at zero bias voltage, it is useful as a benchmark. The 
top left panel of Fig. [2] displays the steady state mag- 
netization predicted by the proposed method at V = 0, 
plotted against the inverse cutoff time jA- at several tem- 
peratures. Comparing the predicted steady state to the 
equilibrium results (circles) shows that the cutoff approx- 
imation converges rapidly even as one crosses the edge 
of the Kondo regime. For the very small magnetic field 
h = 0.0 IT in Fig. [2j the relative errors are rather large, 
but considered on the full scale of the magnetization the 
precision demonstrated here is remarkable. 

The effects of taking the system out of equilibrium are 
illustrated in the lower left panel of Fig. [2] Here a con- 
stant temperature f3T = 1 is maintained while the bias 
voltage V is varied at h = 0.1T; the numerically exact 
V = result is also shown. This plot clearly illustrates 
that convergence of the method generally occurs at even 
shorter times in nonequilibrium conditions — equilibrium 
exhibits the longest memory, consistent with expecta- 
tions. 

An independent approach to verifying convergence re- 
lies on direct examination of individual elements of the 
memory kernel as a function of time. Several representa- 
tive elements are displayed at h = O.Ofr and j3T = 1 in 
the top right panel of Fig. [2] with the inset highlighting 
short times. Within the times accessible by BoldOCA, 
the memory kernel elements go to zero within the nu- 
merical accuracy. Below this, on the same time scale 
and for the same parameters, the time dependence of the 
three distinct elements of the reduced density matrix a is 
plotted for an initially magnetized dot in the lower right 
panel of Fig. [2] With this initial condition and within 
the symmetric Anderson impurity model, the diagonal 
density matrix entries ag and 03, which express charging 
dynamics, are identical. They both relax rapidly, and in 
fact their steady state values could have been obtained 
to very good accuracy within BoldNCA only. The differ- 
ence in scale between the spin relaxation time of o\ , a% 
and the memory decay in the upper panel, however, is 
striking — and is why our memory kernel methods are es- 
sential for obtaining long-time behavior. To obtain a rea- 
sonable converged steady state directly, one would need 
to reach times Tt > 20 with errors of similar magnitude 
compared to what we have obtained at Tt — 2 with the 
current approach. The exponential scaling in time typi- 
cal of all general exact methods makes this unfeasible. 

We now turn from the demonstration of convergence 
to presentation of results. The left panels of Fig. [3] show 
the time evolution of the magnetization from an initially 
polarized state at different voltages and magnetic fields, 
with f3T = 1. At low voltages two separate relaxation 
timescales are apparent: immediate fast relaxation fol- 
lowed by later slow relaxation. At high enough fields 
(bottom), an overshoot effect appears along with oscil- 




Figure 2. Top left panel: The steady-state magnetization 
obtained from the memory formalism at several temperatures 
plotted as a function of the inverse cutoff time, and compared 
in the equilibrium cases to exact CT-QMC results shown as 
circles, for h = 0.01F and V = 0. Bottom left panel: The 
same plot at f3T = 1 and h — 0.1F, for several voltages. Right 
panels: equilibrium memory kernel k (top) and populations <7; 
(bottom) as a function of time for /3T = 1 and h = O.Oir. The 
inset shows the memory kernel at short times. The squares 
in the bottom right panel are approximate OCA results. 



latory behavior which is seen more clearly in the up- 
per right panel. As we increase the voltage the sec- 
ond timescale is suppressed and eventually the relaxation 
becomes exponential. However, the voltage required in 
order to reach this regime is surprisingly large. In the 
top right panel of Fig. [3] we show that nonequilibrium 
NCA and OCA (not supplemented by QMC) do very 
poorly away from h = and cannot be considered reli- 
able, whereas the memory approach and the underlying 
BoldOCA continue to converge very well. 

In the lower right panel of Fig. [3j we show an exam- 
ple of the temperature dependence of the t — > 00 limit 
of the magnetization at constant magnetic field and a 
range of bias voltages. Interestingly, at higher voltages 
(but substantially below ~ w c where the lead chem- 
ical potentials approach the band cutoff) the temper- 
ature dependence becomes non-monotonic. We believe 
this is a population switching effect [33], which leads to a 
suppression of the magnetization by population transfer 
from the magnetized |1) and |2) states to the unmagne- 
tized |0) and |3) states which are activated for V > U. 
The rate for this transfer process is approximately pro- 
portional to the lead occupation at the energy difference 
between the states: / (/?, AE, fi) = 1+e/3( i E _ M) , with AE 
equal to half the interaction energy ^ and /i = y or 
— ¥, depending on the lead. / is therefore an increas- 
ing function of temperature for V < U and a decreasing 
one for V > U . At small voltages the effect of the pop- 
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Figure 3. Left panels: Time dependence of the cutoff- 
converged magnetization at f3T = 1 starting from a fully mag- 
netized dot, at different magnetic fields h and bias voltages 
V. Top right: Comparison with NCA and OCA at h = 2T 
and V — 0. Bottom right: Temperature dependence of the 
h — 0.01F steady state magnetization at different voltages. 
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Figure 4. Generalized magnetic susceptibility \ = T? as a 
function of voltage, for (top) different magnetic fields and 
(bottom) different temperatures, at h = O.ir. Approximate 
results from a master equation calculation are shown in dotted 
lines for the lowest and highest temperatures in the lower 
panel. 



ulation transfer results in a reduction of the magnetiza- 
tion (as expected), while at large voltages the population 
transfer enhances the intermediate-temperature magne- 
tization. At still larger temperatures, the nonequilibrium 
effects are washed away and normal thermal suppression 
of the magnetization occurs. A maximum can therefore 
be expected, as is indeed observed in the numerically 
converged calculations. 

In Fig. [4] we display the steady state voltage depen- 
dence of the generalized magnetic susceptibility x=j[- 
At small h this quantity is /i-independent. The top panel 
shows clearly how the regime in which m is linear in h de- 
pends on voltage at a constant temperature. The bottom 
panel of Fig. [4] shows the voltage dependence at different 
temperatures within the linear regime. One immediately 
noticeable feature is the decrease of \ with increasing /3 
at high voltage, which corresponds to the non-monotonic 
temperature dependence discussed in the bottom panel 
of Fig. |4j A second interesting feature is the fact that 
the plots have a simple, Lorentzian-like structure, sug- 
gesting that the results may be in a regime accessible to 
analytical methods based on performing logarithmic cor- 
rections around rate equations |37j : in the dotted lines 
in the bottom panel we show for comparison results ob- 
tained by solving the classical rate equations (obtained 
by simple perturbation theory to second order in the hy- 
bridization). The large discrepancy between the master 
equation and numerically exact results at j3T — 1 shows 
that in the crossover regime, even at large voltage a cor- 
rect account of the physics requires accurate calculations 
of the sort performed here. 

In conclusion, by unifying numerically exact bold 



Monte Carlo methods with the exact memory approach 
we have developed a new, numerically exact formalism 
free from systematic errors and well suited for the real 
time solution of nonequilibrium quantum impurity mod- 
els. In practice, the capabilities of this formalism are 
unparalleled: the method generates precise, converged 
results at all timescales, in cases where the current state- 
of-the-art approximate methods clearly fail. For the 
nonequilibrium Anderson impurity model, the formalism 
performs well even as one enters the Kondo regime, a 
regime previously inaccessible with accurate numerical 
methods. 

Our formalism has allowed us to explore the detailed 
behavior of the nonequilibrium magnetization, and we 
have made predictions regarding multi-scale, oscillatory 
quenching dynamics at high magnetic fields; the effect of 
voltage on dynamical relaxation; and population-driven 
reversal of the magnetization's temperature dependence 
at high voltages. These results are obtained at param- 
eters where no other currently available method is reli- 
able. As the temperature is further lowered, one expects 
to encounter the formation of Kondo peaks at the chemi- 
cal potential. How this will affect the behavior described 
here remains an interesting and open question, and work 
is currently being carried out to further investigate this 
issue. Future research will address lower temperatures 
and a wider variety of observables; it is also worth stress- 
ing that both bold techniques and the memory formalism 
are not specific to the Anderson impurity model, and are 
expected to have many more applications. 
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